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Abstract. The diagonal case of the sl(2) Richardson-Gaudin quantum pairing model 
^ El El EI El El d is known to be solvable as an Abel-Jacobi inversion problem 
[lISinBHHEllIHimilHl. This is an isospectral (stationary) solution to a more general 
integrable hierarchy, in which the full time evolution can be written as isomonodromic 
deformations. Physically, the more general solution is appropriate when the single- 
particle electronic spectrum is subject to external perturbations. The asymptotic 
behavior of the nonlinear oscillations in the case of elliptic solutions is derived. 



1. Introduction 

The integrable system described by the pairing hamiltonian introduced by Richardson 
and Sherman [H El El IH El Ej in the context of nuclear physics has received revived 
interest in recent years, after being applied to metallic superconducting grains (THj and 
cold fermionic systems [13 HH] . The model is intimately related [TBI HH1 E0j to a class 
of integrable systems generally referred to as Gaudin magnets [Jj. These systems have 
been studied both at quantum and classical level [HlESEHlEIlEElEIIlEZlEElEniEn!, 
in the elliptic well as trigonometric and rational degenerations, using various 

methods from integrable vertex models to singular limits of Chern-Simons theory. 

In an interesting regime of the pairing problem was considered, which may be 
relevant to recent experiments with cold fermionic gases exhibiting the paired BCS state 
[SU EH EH EH EE1 EH E3 EHj • It was shown that for such systems, the time scales of the 
order parameter ta ~ | Z\ | 1 , and the quasiparticle energy relaxation time r e are both 
much larger than typical time for switching on the pairing interaction r , essentially 
given by the variation of external parameters, such as detuning from the Feshbach 
resonance. It was argued that in this regime, for times t C r £) the dynamics of the 
system is given by non-linear, non-dissipative equations describing the coherent BCS 
fluctuations for the system out of equilibrium. In this limit, the system is integrable, 
and features non-perturbative behavior, such as soliton-type solutions. 
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In the mean-field limit, such non-trivial solutions describing the collective mode of 
the Anderson spins were derived in J7| , for a two-level effective system. This work 
was generalized ^Hj in algebro-geometrical terms. 

In (401 El H2j; the long-time behavior of the solution has been considered, under 
various conditions. An issue not addressed so far is the relaxation of the nonlinear 
oscillatory solution induced by perturbations of the spectral curve, physically justified by 
coupling to the environment. Several possible kinds of perturbations may be considered, 
which may lead to different types of relaxation. 

In this paper, we consider the effect of fluctuations of single-particle energy levels, 
which amount to slow deformations of the Liouville tori, and can be described by 
hydrodynamic-type equations in phase space O Ej. These equations describe the 
evolution of the moduli for the complex curve of the system. 

2. The Richardson-Gaudin Model 

2.1. The quantum pairing hamiltonian 

Following [21], we briefly review the Richardson pairing model. It describes a system 
of n fermions characterized by a set of independent one-particle states of energies q, 
where the label I takes values from a set A. The labels may refer, for instance, to 
orbital angular momentum eigenstates. Each state I has a total degeneracy di, and the 
states within the subspace corresponding to I are further labeled by an internal quantum 
number s. For instance, if the quantum number I labels orbital momentum eigenstates, 
then di = 21 + 1 and s = —I, . . . ,1. However, the internal degrees of freedom can be 
defined independently of I. We will assume that di is even, so for every state (Is), there is 
another one related by time reversal symmetry (Is). For simplicity, we specialize to the 
case d\ — 2, s =|, |. Let c\ s represent the fermionic creation operator for the state (Is). 
Using the Anderson pseudo-spin operators (HH) (quadratic pairing operators), satisfying 
the su(2) algebra 

[tl tf] = ±%t± [t+ tj] = 25 i5 t), (1) 
the Richardson pairing hamiltonian is given by 

Hp = 2 ^f -gJ2 = E 2eit * - 9 t+ • t_ ' ( 2 ) 

zeA i,i' zeA 

where t = Yli tz is the total spin operator. It maps to the reduced BCS model 

H = Yl £ p a k- 6 P^ - 9 Yl cp T C-p|C-k|C kT (3) 

p,cr p,k 

by replacing the translational degrees of freedom by rotational ones, where I G A = 
{l,...n} ennumerates the one-particle orbital degrees of freedom, while s =|, j 
indicates the two internal spin states per orbital (di = 2). The pairing hamiltonian 
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can be decomposed into the linear combination 



H p = 2J2^iRi + 9 



leA 



f3 



(4) 



At a fixed value of the component t of the total angular momentum, the last term 
becomes a constant and is dropped from the hamiltonian. The operators Ri (generalized 
Gaudin magnets jjj) are given by 



Ri 



t 3 - 



g ^ t r t v 



E 



2 14 — ' 6; — £;/ 



(5) 



These operators solve the Richardson pairing hamiltonian because ^B] they are 
independent, commute with each other, and span all the degrees of freedom of the 
system. Richardson showed PUE] that the exact N— pair wavefunction of his hamiltonian 
is given by application of operators b\ = ^ 2e l ~e k ^° vacuum (zero pairs state). The 
unnormalized N— pair wavefunction reads ^(e,) = Ylk=i ^ll^)- The eigenvalues e k 
satisfy the self-consistent algebraic equations 



1 

9 



1 



2e z - e k 



which can be given a 2D electrostatic interpretation [21] with energy 



U(e h e k ) 



2 
9 



N 



Lfe=l 



1=1 



+ 



(6) 



(7) 



A? 



2^^1og|e fc - 2ei\ - 4 ^ log |e fc - e p \ - lo S l 2e ' ~ 2e il ( 8 ) 

i=l fc=l fc<p i<j 

Equations © appear as equilibrium conditions for a set of charges of strength q = 2 
placed at points e k , in the presence of fixed charges of strength q — — 1 at points 
2e/, and uniform electric field of strength -, pointing along the real axis. This 
interpretation proves to be very useful for the conformal field theory (CFT) description 
of the Richardson problem. The electrostatic energy (JJJ) is minimized for values {e k } 
corresponding to pair energies. In (|7|). n, N represent the number of single-particle 
levels and the number of pair energies, respectively. For large interaction constant g, 
the equilibrium positions {e k } form a set of complex conjugated pairs defining a curve 7 
in the complex plane of energies. We note that the eigenvalues of Gaudin hamiltonia 
are proportional in this language to the values of the electric field at positions 2e.- l) 

For a set of single-particle energies {e^}, the BCS ground state is obtained by 
minimizing the electrostatic energy (j7J) with respect to positions of the free charges at 
{e k }. Once found, they also determine exactly the values of the electric field at positions 
{2ej} on the real axis, which are proportional to the ground-state eigenvalues {rf 5 }. For 
any other values of {r^}, the electrostatic energy (JJJ) is not minimized. This indicates 
that for arbitrary values 7^ rf s , the system is not in equilibrium. 
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2.2. The mean field limit of Richardson- Gaudin models 

General description of the classical model In the mean-field limit, the spin operators 
ti are replaced by their quantum mechanical averages. Written in terms of the classical 
vectors S; = 2(t;), the semiclassical approximation for the pairing hamiltonian becomes 

iW = J>S?-||J-| 2 , (9) 

zeA 

where J = J2ieA^i aric ^ the BCS gap function is given by A = gJ~/2. Replacing 
commutators by canonical Poisson brackets, 

{S?, Sf} = 2e a ^S78 ij , (10) 

variables Sf become smooth functions of time. In this limit, the problem can analyzed 
with tools of classical integrable systems, and the solution is known to be exact as 
n — > oo. 

The Poisson brakets (|10|1 and hamiltonian Q lead to the equations of motion 

Si = 2(-A + e. t z)x S h (11) 

where 2A = (gj x ,gj y ,0) and J is the total spin. The semiclassical limit of Gaudin 
hamiltonia are independent constants of motion, 

h = 0. (12) 



qz 9 \^ Si ■ Sj 
^ ~ 2 2_, 7~^7- 



Equations (fTTj) describe a set of strongly interacting spins and have generic non-linear 
oscillatory solutions. The exact solution may be obtained through the Abel-Jacobi 
inverse map [01 ESI HHj • 

This solution can be described exactly in the language of hyperelliptic Riemann 
surfaces (see Appendix A for details). At this point, it is useful to make use of the 
intuitively clear features of this construction (Figure Q). For a given set of initial 
conditions for the spins {Si}, i.e. also of the constants of motion {r^}, a polynomial 
Q{u) of degree 2n and with n pairs of complex conjugated roots E 2 k+ 2 = E 2 k+i, k = 
0, . . . , n— 1, is constructed. A schematic representation of these roots is given in Figure^ 
Between each pair of roots, we place a simple cut Ck = [E 2 k+i, E 2 k+2\ on the complex 
plane of energies. The surface thus obtained is a representation of a torus of smooth 
genus g — n — 1. 

Variables are introduced for n — 1 of these cuts, with respect to which the 
equations of motion separate. The variables evolve in time in a complicated fashion, 
solving a system of nonlinear coupled differential equations ()A.4jl . Up to a constant, the 
time dependence of the gap parameter amplitude is given by 

log|A(t)|=Tm u(t)dt, u{t) = ^u k {t). (13) 

** k 

The widths of the cuts Ck and the periods of the non-linear oscillators Uk are determined 
by the values of constants of motion {r^}. For the particular choice = rf , all 
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Figure 1. Schematic representation of the Liouvillc torus for the integrable system 
CD- 



the cuts Cfc, k = 1, ... ,n — 1 vanish, and the width of the remaining cut C n equals 
the equilibrium value of the gap function: {E^f-i — E^l — 2|A| G5 . In that case, 
the oscillators Uf. = E 2 k-i = E 2 k are at rest, and the only time dependence left in 
the system is the uniform precession of the parallel planar spins S~ , with frequency 
uj = 2 Ylk=i e k ~~ 2 Ylp=l Ei v -\ — Ys7=i ^i- m the case of particle-hole symmetry, u 
vanishes as well. 

2.3. The onset of pairing interaction and long-time behavior of oscillations 

Equations do not impose any particular constraints on the integrals of motion. 
These constants are known for the metallic state t < 0, but change abruptly during 
the short interval when the electron-electron interaction is turned on. The onset of 
pairing interaction is a delicate problem in itself |%U1 HT| H2*] . which deserves further 
investigation. There are several reasons which make the issue non-trivial. We review 
here several of them. 

a) In the case of KdV and KP2 hierarchies, it is known that finite-gap solutions (for 
which only a few cuts Ck are non-degenerate) are dense in the space of all periodic 
solutions. A conjecture of Krichever |32] extends this fact to arbitrary algebraic 
curves. Therefore, solutions where an infinite number of independent frequencies 
contribute to the mean-field approximation require very special initial conditions. 
The system is most likely to be described either by (i) a finite-gap solution, or 
(ii) a non-oscillatory function. The case (ii) has not been investigated in previous 
studies. 

b) Another issue related to solutions described by only a few frequencies is that the 
quantum integrability of the system may impose additional constraints and provide 
a selection criterion, non-existent in the mean-field approximation. The exact 
solution of quantum XXX Gaudin system under a sudden variation of the spectral 
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curve is a subject of active research and will be addressed elsewhere |%o] . 

c) Concerning the finite-gap solutions of the mean-field approximation, we note that 
the onset of attractive interaction is realized through some mechanism which 
effectively changes the parameters of the spectral curve. It is reasonable to 
assume that this variation is not instantaneous and in fact may continue as a 
perturbation for the rest of the evolution of the system. This mechanism describes 
coupling of the system to the environment (the external magnetic field and optical 
trap, for degenerate Fermi gases) and therefore may induce fluctuations in the 
parameters describing the spectral curve. We therefore consider the effect of such 
perturbations which do not change the monodromy of the solution. The concept 
of isomonodromic deformations of integrable systems was introduced in |4*T1 |4*H1 H§] 
and used extensively in the theory of Painleve transcendents. We give here a brief 
review of the method. 

We consider the unperturbed problem given by a nonlinear scalar differential 
equation R(u,u,u, . . .) = 0, where R is a rational function of the solution u(t) and 
its derivatives u,u,.... Specializing to the cases described by hyperelliptic spectral 
curves, we express it through the Lax pair of 2 x 2 operators L, A, solving the linear 
vector problem 

L$ = (M, d t ^ = AV, d t L = [A, L], (14) 

where A,/x are two auxiliary complex variables, A(X,u) a 2 x 2 matrix, and \I/(A,/i, t) is 
a column vector. The operators are chosen such that elements of the identity L = [A, L] 
are equivalent with R(u, ii, ii, . . .) = 0. A trivial calculation gives the isospectral property 
fi = 0. 

Specializing the operator L to have the form [^0] (the u functional dependence is 
implicit throughout) 

™=[S <*> 

the eigenvalue equation for \1/ becomes 

fi 2 + det L = 0, fi{X) = ±VdetL(A). (16) 

When functions a, b, c in (JT5|l are rational, equation (fTB^l defines the hyperelliptic 
Riemann surface /i(A) called spectral curve of the system ()14j) . As noted before, the 
time evolution leaves the spectral curve invariant, which is why the problem is sometimes 
called isospectral. 

The isomonodromic deformation [H] is introduced through 

m = ^ + ed x V, [ed x -L,dt-A]= 0, (17) 

where the second equation is the isomonodromy requirement, and 0<e<l. We recover 
the unperturbed problem by setting e = 0. The first equation is easily integrated and 
gives the formal solutions 



log*±(A,/i,t) = - 



-//A ± % J a/ det L{a)da 



(18) 
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up to constants in A. Let us now impose the saddle-point (or turning-point) condition 
d\\og^> = 0. This will give /x = det L, i.e. the spectral curve. Therefore, we 
may see the isospectral problem as a saddle-point (turning point) approximation. The 
compatibility equations 

[ed x - L, d t - A] = (19) 

can be recast in the form 

e(d v L - d\A) + d T L - [A, L] = 0, (20) 

where we have split the time dependence d t into a fast time scale d T and a slow one ed v . 
At zero order in e, (|2())l is simply the unperturbed problem. The first order correction 
gives the slow-time scale dependence of the modulated solution u e (r, v). Moreover, from 
the matrix elements of (12 OH . we get after averaging over the fast motion in r |51j . 

d v detL = -(2aA 11 + bA 21 +cA 12 ), (21) 

where the bar signifies r— averaging. Since for the unperturbed problem, d T det L = 0, 
the averaging is justified. 

Equation (}2T|) is simply the result of Bogoliubov-Whitham [52] averaging for the 
problem ()17j) . It tells us how the previously invariant spectral curve now changes slowly 
in time over the large time scale v. It also gives us the form of the deformed solution 
u e (r, ki(v)), as a modulation of the original solution. 

Note The fact that such perturbations generalize the autonomous Gamier system of 
[To] to a non-autonomous system of Schlesinger type was indicated in j^JJ. Extensions 
of these systems to include a constant matrix were discussed in jS3J. . 

3. The effect of weak perturbations 

3.1. Topological classification of finite-gap solutions 

In the presence of spectrum symmetry e k = — e_fc, SI = —S^_ k , the distribution of cuts 
Cfc obeys the same symmetry. The number of non-degenerate cuts is therefore even. 
This analysis uses the fact that at t — all Anderson pseudo-spins are aligned along 
the z axis. In the simplest non-trivial case, there are only two non-degenerate cuts as 
shown in Figure El and g = 1, while the corresponding variable u\ is given by elliptic 
functions. Other non-trivial cuts Ck may exist in general, associated with dynamics of 
variables Uk- In the limit of "small" cuts [SI], their contributions separate are given by 
trigonometric functions and the behavior of the gap parameter takes the simplified form 

iog jsjH =im I u[t)dt = Ueu[t) + ^ (22) 

As we shall see, the most relevant contribution is due to the elliptic part U e u, studied 
in the next section. 
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Figure 2. Distribution of complex roots Ei around the origin, for the elliptic 
approximation. The mean level spacing is d. 



3.2. Modulations of elliptic solutions 

For the root distribution shown in Figure there is only one variable Uk, taking 
imaginary values [TH]. It solves an equation of the type [TH] 

(ii) 2 + (u 2 + m 2 )(u 2 + M 2 ) = 0. (23) 

For simplicity, we therefore make the transformation u — > iu, keeping the time variable 
real. The new function is a real quartic oscillator which satisfies 

(u) 2 = (u 2 - M 2 )(u 2 -m 2 ). (24) 

The solutions corresponding to this distribution of roots are 

u 1 (t) = m-sn(Mt + <f> u k 1 ), k x = m/M, (25) 
u 2 (t) = M ■ sn(mt + (f) 2 ,k 2 ), k 2 = M/m, (26) 

where sn is the Jacobi sine function, and <pi t2 are arbitrary phases. In the degenerate 
case m = M, the solutions become hyperbolic functions. Solution u 2 is non-physical in 
our case. 

In order to set up the isomonodromic deformation method, consider the Lax pair 

L = — [2uax + Au\a 2 + (4A 2 - £ + 2u 2 )ia 3 ] , (27) 
A = i\a 3 - ua 2 , ed x ^ = LV, d t V = A®, (28) 

where o~ a , a = 1, 2, 3 are the Pauli matrices, e is a small real number, and \& = (^i, V^)* 
is the Baker- Akhiezer function. Choosing the initial condition ipx = 0, ip 2 = I, we can 
extract the amplitude of the ratio A(t)/A(0) from ty(t) as 

A(t) [Mt) + Mt\ (29) 



A(0) 2 
The compatibility (zero-curvature) conditions 

[d t - A, ed x - L] = (30) 
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dti = e, d 2 t u = 2w 3 - fit. 



We note that 



detL 



-4{(n) 



^ 2 + (2A 2 -£/2) 2 ]}. 



(31) 



(32) 



Setting e = gives the equation (u) 2 — u A + (,u 2 = constant, while the limit e = 1 yields 
the Painleve II equation. In fact, the unperturbed case e = allows to retrieve the full 
elliptic solution, from the equation 

= 0, det(L) = 0, (33) 

which gives the elliptical function u satisfying 



u 



[u 4 - £u 2 + (2A 2 - e/2) 1 



0. 



(34) 



The physical solution u\ is obtained identifying 

£ = m 2 + M 2 , 2A = M + m. (35) 

Now restore e ^ 0, write d t = d T + e<9g, and retain terms of order e from f)3Uj) . Averaging 
over the fast variable r gives 

det L = L 2 2d x A u + LudxA2 2 , (36) 

or equivalently, 

d^detL = -(4u2- 2^ + 8A 2 ). (37) 



Performing the computations, we obtain 

d^-^u 2 - (u) 2 } = -1?. 

Writing the elliptic equation as 

[uf = u 4 - £u 2 + /i 2 , 
the physical solution takes the form 



(38) 
(39) 



ui(t) 



sn 



up to an arbitrary initial phase 0, and elliptic modulus 
1 - 



k 



1 



:f) 2 



i + 



i-(f) 2 



The Whitham averaging equation has the form 
,<V 1 2 f(P) 



4- 



1 



(40) 



(41) 



<9£ 2 2-k 2 1C{k 2 y 
where £, /C are the complete elliptic integrals of the first and second kind, respectively. 
Equation (|41|) has a fixed point at k = 0, | — > 0. This shows that, on the slow 
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time scale, the parameter ///£ = m/M + 0(m 3 /M 3 ) goes to zero, as £ = m 2 + M 2 
increases. Expanding the solution in this limit, and integrating under the separation of 
time scales asumption, we obtain for the elliptic contribution to the gap parameter, the 
approximation 

A eU (t) = A(0)e k f sn(T > k2)dT , (42) 
and k 2 = {£ + 0(/i 4 /«e 4 ) -> as t -> oo. 

Asymptotic behavior of modulated elliptic solutions Starting from the Lax pair (|27| 
ES|I . we can obtain the asymptotic behavior of \I/(A,t) as t —>■ oo, in the whole complex 
A plane [55J. The analysis is simpler when working in the variable z = There are six 
Stokes sectors at z — > oo, with canonical asymptotes for but the region of interest is 
z — 1/2, where 

4A 2 /x 2 , x 

T -l. (43, 

Using the Whitham equation, separation of scales, and identifying the strength of the 
fluctuations e with an effective temperature T, we obtain for the elliptic contribution at 
large times (J22J) 

Ueu(t) = (44) 
where T is a bounded oscillatory function \T\ = 0(1). 

Modulations of trigonometric solutions In order to analyze the slow dynamics of the 
small cuts, let £ — > oo in an d write the solution as 

£. 

This assumption is consistent with Formally, this is simply the trigonometric 

approximation of the general elliptic solution, in the limit of small modulus k 2 . However, 
the parameter £ in (|45jl must be sent to oo mush faster than the physical £ = M 2 +m 2 , in 
order to give the correct spherical limit. The contributions from the small cuts therefore 
vanish faster than the elliptic component. 



// =-9cos( v / £t + 0), I^V^O. (45) 



Higher-genus contributions In the case where there are several non-degenerate cuts at 
t — corresponding to a high genus g > 1, the isomodromic deformation method can 
be applied in the same way as for g — 1, leading to modulated hyper-elliptic solutions of 
Painleve equations degree higher than 2. There are few systematic results in this field, 
and a complete classification does not exist at this time. Specific high-genus solutions 
may have interesting topological properties with physical interpretations in terms of the 
collective excitations of the Gaudin magnet |46j . 



Relaxation of nonlinear oscillations in BCS superconductivity 



11 



Appendix A. Isospectral case and the Abel-Jacobi inversion problem 



In |3 , the system (fTTjl with fixed spectral curve was solved through inverse Abel- 
Jacobi mapping, by using Sklyanin separation of variables techniques [El EE]- Interesting 
connections to generalized Neumann systems and Hitchin systems were discovered in 
[21]. The solution starts from the Lax operator 



^ / , n 2 v — \ Sa ■ a 



i=0 



a(X) 6(A) 
c(A) -a(A) 



(A.l) 



where a a , a = 1,2,3 are the Pauli matrices, and A is an additional complex variable, 
the spectral parameter. Let Uk, k = 1, ... ,n — 1 be the roots of the coefficient c(A). 
Poisson brackets for variables 5" read 

{Sf, Sf } = 2e a ^S25 jk (A.2) 

The Lax operator (|A.1|) defines a Riemann surface (the spectral curve) T(y, A) of genus 
g = n — 1, through 

p(\y 2 



Q(A) = det£(A) 



where P(A) = n^(A - <*). 

The equations of motion for the hamiltonian Q become 



(A.3) 



U; 



2iy(ui 



m-1 



iJ- 



J- 



U; 



U ; 



gJ 3 + 2j2tk 



u 



fe=i 



(A.4) 



In (Q, u = -2 M ^ K^) = o. 

From the equations of motion, it is clear that knowledge of the initial amplitude 
of J~ and of the roots {ui} is enough to specify the n unit vectors {Si}, for a given 
set of constants of motion {Ri} given by the classical limit of Gaudin hamiltonia. The 
Dubrovin equations (jA.4|) are solved by the inverse of the Abel-Jacobi map, as we explain 
in the following. We begin by noting that the polynomial Q(X) has degree 2n, and is 
positively defined on the real A axis. Therefore, the curve T(y, A) has n cuts between 
the pairs of complex roots [E 2 i-i, E 2 i),i = 1,2, ... ,n, perpendicular to the real A axis. 
The points Ui belong to n — 1 of these cuts, Ui G [E 2 i-i, E 2 %\,i = 1, ... ,n — 1. These 
g = n — 1 cuts allow to define a canonical homology basis of T, consisting of cycles 
{ai,j3i},i = \,...,g. With respect to these cycles, a basis of normalized holomorphic 
differentials {uJi} can be defined, through 



y 



Mi, 



Hi 



M 



-1 



(A.5) 



The period matrix Bij = Jg OUi is symmetric and has positively defined imaginary part. 
The Riemann 9 function is defined with the help of the period matrix as 



6(z\B) = e 



2m(n t z+^n t Bn) 



(A.6) 



nez9 
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The g vectors Bj, consisting of columns of B and the basic vectors define a lattice in 
C 9 . The Jacobian variety of the curve T, is then the g— dimensional torus defined as the 
quotien J(T) = C 9 / (Z 9 + BX 9 ). The Abel-Jacobi map associates to any point P on T, a 
point (g— dimensional complex vector) on the Jacobian variety, through A(P) = f^<*>- 
Considering now a g— dimensional complex vector of points {Pk},k = l,...,g on T, 
defined up to a permutation, we can associate to it the point on the Jacobian 

9 

z = a(P) = Y,A(Pk) + K, (A.7) 

fc=i 

where K is the Riemann characteristic vector for T. 

The map (jA.7|) suggests that we now have a way to describe the dynamics on T 
by following the image point on the Jacobian. Given a point on the g— dimensional 
Jacobian z = (d, . . . , Cn-i), we can find an unique set of points {A^}, k = 1, . . . , g on V, 
such that z = a(A), and 9(a(P) — z\B) = 0. The system evolves in time according to 
the point z(t) 

(k = ic k , 1 < i < 9 ~ 1, Cn-i = i(°n-i + t), (A. 8) 

where {c^} is a set of initial conditions, such that zq = zit = 0) = a(c), and c is the set 
of initial conditions for positions of A on T. Together with the initial condition which 
determines the initial amplitude of J~ , this set will determine entirely the evolution of 
the functions Ui(t), J~(t). 
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